Modulational instability of Rossby and drift waves and generation of zonal jets 
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We study the modulational instability of geophysical Rossby and plasma drift waves within the 
Charney-Hasegawa-Mima (CHM) model both theoretically, using truncated (four-mode and three- 
mode) models, and numerically, using direct simulations of CHM equation in the Fourier space. The 
hnear theory [Ij predicts instabiUty for any amplitude of the primary wave. For strong primary waves 
the most unstable modes are perpendicular to the primary wave, which correspond to generation 
of a zonal flow if the primary wave is purely meridional. For weaker waves, the maximum growth 
occurs for off-zonal inclined modulations. For very weak primary waves the unstable waves are 
close to being in three-wave resonance with the primary wave. The nonlinear theory [5] predicts 
that the zonal flows generated by the linear instability experience pinching into narrow zonal jets. 
Our numerical simulations confirm the theoretical predictions of the linear theory as well as of the 
nonlinear pinching. We find that, for strong primary waves, these narrow zonal jets further roll up 
into Karman-like vortex streets. On the other hand, for weak primary waves, the growth of the 
unstable mode reverses and the system oscillates between a dominant jet and a dominate primary 
wave. The 2D vortex streets appear to be more stable than purely ID zonal jets, and their zonal- 
averaged speed can reach amplitudes much stronger than is allowed by the Rayleigh-Kuo instability 
criterion for the ID case. We find that the truncation models work well for both the linear stage and 
and often even for the medium-term nonlinear behavior. In the long term, the system transitions to 
turbulence helped by the vortex-pairing instability (for strong waves) and by the resonant wave- wave 
interactions (for weak waves). 



PACS numbers: PACS go liere 

INTRODUCTION AND MOTIVATION 

Zonal flows are prominent features in the atmospheres 
of giant planets such as Jupiter and Saturn [H HI |5] , the 
Earth's atmosphere [B] and its oceans [3 |7]. Geophys- 
ical jets can regulate small-scale turbulence and trans- 
port processes via, for example, the "Barotropic Gover- 
nor" mechanism [S]. Zonal flows are also important in 
plasma turbulence of fusion devices p^ . There, they can 
also regulate the turbulence and suppress transport via a 
drift- wave / zonal-flow feedback loop [101 CI] ■ The latter 
process is presently considered the main mechanism for 
the Low-to-High (LH) confinement transitions in toka- 
maks discovered in |12j . - an effect which is crucial for 
the success of future fusion devices. 

Two main zonal flow generation scenarios have been 
discussed in the literature. According to the first scenario 
zonal flows are generated via an inverse energy cascade, 
which could be local or nonlocal [TT] . The mechanism 
for such an inverse cascade is similar to that of 2D Navier- 
Stokes turbulence [13], but the presence of the beta-effect 
makes this cascade anisotropic. This leads to a prefer- 
ential transfer of energy into zonal flows at large scales 
rather than into round vortices as would be the case in 
Navier-Stokes turbulence. The beta-effect leads to three- 
wave resonant interactions which preserve an additional 
(to the energy and the potential enstrophy) quadratic 
invariant, - zonostrophy [TU [TSJ US]. Application of the 



standard Fj0rtoft argument to the three invariants, the 
energy, potential enstrophy and zonostrophy lead to the 
conclusion that the energy can only be transferred to 
large zonal scales [13]. This statistical argument is ex- 
plained in detail in |17j . The second mechanism of zonal 
flow generation, and the principle topic of this article, 
is via modulational instability of a primary meridional 
Rossby/drift wave. In practice, such primary waves are 
themselves the result of an instability (typically the baro- 
clinic instability in GFD or the ion-temperature-gradient 
instabihty in tokamaks) [TJ d HHJ (H UHl HIl 112 ■ 

These mechanisms are unlikely to be exclusive in prac- 
tice and both may coexist under some conditions. The 
extent to which one mechanism dominates over the other 
is determined by the parameter regime and configu- 
rational details. If the parameter regime were to be 
such that the baroclinic instability resulted in meridional 
Rossby waves, zonal flows would presumably result from 
the MI mechanism, whereas if the parameter regime were 
to be such that the baroclinic instability resulted in more 
isotropic eddies at the Rossby deformation radius, the 
cascade scenario would likely be more relevant. In our 
purely barotropic model, these effects are modeled by the 
initial condition. A narrow initial spectrum of the waves 
and large initial amplitude promotes the modulational in- 
stability mechanism leading to fast zonal flow generation 
bypassing the turbulent cascade stages. On the other 
hand, for broad initial spectra, the cascade scenario is 
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likely to be more relevant. There is an analogy with the 
turbulence of surface gravity waves on water where the 
inverse cascade and modulational (Benjamin-Fair) insta- 
bility [23| can compete with each other in the generation 
of long waves [53]. A quantitative measure, called the 
Benjamin- Feir index, was suggested to estimate proba- 
bility for triggering the modulational instability [2H I25j . 
Developing a similar approach for the Rossby/drift wave 
system would also be useful. However, we will leave this 
interesting subject for future studies, and in the present 
paper we will only be concerned with the modulational 
instability of a monochromatic wave. 

We will start by revisiting the linear theory of the mod- 
ulational instability which was first analysed by Loretz 
[18] and then treated in great detail in a beautiful pa- 
per by Gill [T]. Using numerical and semi-analytical cal- 
culations we highlight the most important properties of 
Gill's theory. In particular, we will see how the charac- 
ter of instability changes with the strength of the carrier 
wave: from being the classical hydrodynamic instability 
of the sinusoidal (Kolmogorov) shear flow for large am- 
plitudes to becoming a (three-wave) decay instability 
of weakly nonlinear waves for small amplitudes [27j . We 
will also study the effect of the finite Rossby/Larmor ra- 
dius on the instability. 

We will then proceed to study the nonlinear stage of 
the modulational instability with direct numerical simu- 
lations (DNS), comparing them with the solutions of the 
four-mode truncated (4MT) and the three-mode trun- 
cated (3MT) systems. We find that at the nonlin- 
ear stage, for strong primary waves, the growth of the 
zonal mode deviates from the truncated dynamics and 
the zonal fiow tends to focus into narrow jets, as was 
theoretically predicted in [5]. These zonal jets subse- 
quently become unstable and acquire the interesting two- 
dimensional structure of a double (Karman-like) vortex 
street. The vortex street appears to be more stable than 
a plane parallel shear fiow with the same zonal profile 
[2H] and persists for a relatively long time until (possi- 
bly due to dissipation) a vortex pairing instability sets in 
and triggers a transition to turbulence [2S]. As the non- 
linearity of the primary wave is decreased we find that 
there is a level of nonlinearity below which this sequence 
of events changes. For sufficiently weak primary waves, 
the jet strength reaches a maximum which is still stable. 
After this maximum is reached, the jet amplitude starts 
decreasing again, continuing to follow the truncated dy- 
namics, and avoids the roll- up into vortices. This re- 
versal of the jet growth, particularly the maximum jet 
strength, is well predicted by nonlinear oscillatory solu- 
tions of the 4MT, and often by the 3MT, equations. The 
latter are relevant for non-degenerate (in a sense which 
we shall explain) resonant wave triads. Once the full sys- 
tem deviates from the solutions of the truncated system, 
as it inevitably does, it sometimes continues to exhibit 
oscillatory behaviour for a while in the weak nonlinear- 



ity cases. These subsequent oscillations have different 
periods, however, and are often rather irregular. 

Along the way, we will examine the relative perfor- 
mance of the 3MT vs 4MT models thereby clarifying 
possible confusions on whether the principal mechanism 
of the modulational instability is 3- wave or 4- wave. 



THE MODEL 

Geophysical and plasma zonal fiows are often men- 
tioned together because of the same simplified nonlinear 
PDE which was suggested for their description, namely, 
the Charney-Hasegawa-Mima (CHM) equation [^ [50] : 

9t(A^-FV;)-H/39.V^ + J[V',AV'] = i^„(-A)"V, (1) 

where "0 is the streamfunction, _F = 1 /p^ with p be- 
ing the deformation radius in the GFD context and the 
ion Larmour radius in the plasma context, /3 is a con- 
stant proportional to the gradient of the horizontal rota- 
tion frequency or of the plasma density in the GDF and 
plasma contexts respectively. We introduced notation for 
the Jacobean operator. 



J[f.9] = {dJ){dyg)~{dyf){d.g). 



(2) 



In the GFD context, the x-axis is in the west-east and the 
y-axis is along the south-north directions respectively. In 
plasmas, the y-axis is along the plasma density gradient 
and the x-axis is, of course, transverse to this direction. 
Also, keeping in mind our numerical simulations which 
will be described below, we introduced to the right-hand 
side (RHS) a hyperviscous dissipation of some degree n > 
2 (a positive integer) and a small positive coefficient Vn- 
Introducing the Fourier transform of the streamfunc- 
tion, V'k = / V'(x)e-*('^-'') dx, the CHM equation, Eq. ([l]), 
ignoring the hyper-viscosity term for now, is equivalent 
to the following: 

+ J ^ T(k,ki,k2)^kiV'k.5(k-ki+k2), (3) 



ki,k2 



where 



T(k,ki,k2) = 



(ki X k2), {kl 



k^ + F 



(4) 



(5) 



and k = {kx,ky) and k = |k|. The system clearly sup- 
ports linear waves, known as Rossby or drift waves, in 
the GFD and plasma contexts respectively. They have 
the anisotropic dispersion relation given by Eq. The 
structure of the nonlinear interaction, Eq. ([5]), is such 
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that the nonlinear term vanishes for a monochromatic 
wave. Hence Rossby waves are actually exact solutions 
of the full CHM equation. Much of this article will focus 
on the stability properties of these solutions. 

Originally, the waves arise due to a primary instabil- 
ity, e.g. the baroclinic instability in GFD PH] or the ion- 
temperature-gradient instability (ITG) in fusion plasmas 
[31]. The instability is not included in the CHM equa- 
tion, and it could be modeled by simulating CHM with an 
initial condition or introducing a linear forcing term on 
the RHS mimicking the linear instability (this would not 
take into account the nonlinear mechanisms in the wave 
forcing). It is interesting that the GFD-plasma analogy 
extends to the instabilities too in that the most unstable 
mode is "meridional" (i.e. along the x-axis) and con- 
centrated at the scales of the order of p. Thus, in most 
of our considerations below we will consider the initial 
(primary) wave which is purely meridional. 

SPECTRAL TRUNCATIONS 

We shall use spectral trunctions of Eq. ([3| in our study 
of the stability properties of Rossby waves. They provide 
approximations of an intermediate degree of complexity 
between monochromatic waves and the full PDE. At this 
stage, such truncations should be viewed as ad-hoc since, 
in reality, all triads are coupled together in Eq. ([s]) . Their 
usefulness will determined by comparision with DNS so- 
lutions of the full system, Eq. ([3]). We shall consider 
two natural truncations: the 3-mode truncation and the 
4-mode truncation. 

3-Mode Truncation (3MT) 

The simplest such truncation is to restrict the RHS of 
Eq. (§ to a single triad containing only 3 modes which 
we shall denote by p, q and p_ = p — q. We construct 
the truncated equations be taking each wave vector in the 
triad in turn and assigning it to be k in Eq. ([s]), enumer- 
ating all ways of assigning the others and their negatives 
to ki and k2 on the RHS and neglect all terms which 
involve "^k's within the triad. Since V'k is the Fourier 
transform of a real field, -^-k — '0k- We then arrive at 
the following equations for the 3-mode truncation: 

dt^p + iLOpijjp = T(p, q,p_)'i/'q^p_ 

atV'q + it^qV'q ^ ^(q, p, -p_ )l/'pV'p_ (6) 

i9t?Ap_ Wp_ '0p_ = T(p_,p, -q)i/'p-0q. 

For most of what follows, it will be convenient to deal 
with Eqs. ([6]) in the interaction representation. Introduce 
^'i^(t) = ?/;k(t)e"''^''*. In terms of ^-k, Eqs. ^ become 

dt^p = r(p,q,p_)*q«'p_e^^-* 



at^Tq = T(q,p,-p_)vl'pvl-p_e-'^-* (7) 
5t*p_ = T(p_,p,-q)vl'p^qe-''^-*, 

where 

A_ = CJp - - Wp_ . 

A similar set of equations can be derived for the other 
natural triad, (p, — q, p+) where p+ = p -I- q: 

5t*p = r(p,-q,p+)^q^'p_e'^+* 
5**q = T(q,-p,p+)^pvl'p^e^^+* (8) 
dt-^^^ = r(p+,p,q)vI/pVl/qe-^^+*, 

where 

A+ = CJp + tJq - LOp^ . 

If A+ = 0, the triad is exactly resonant. Then Eqs. ([s]) 
form an exactly integrable set of equations which have 
been extensively studied [32l [33] . 



4-Mode Truncation (4MT) 

The 4MT model retains both triads, (p,q, p+) and 
(p, — q, p_), where p± = p±q mentioned above. A 
thorough analysis of the 4-mode truncation for the Gen- 
eralised Hasegawa-Mima equation in the case of weak 
nonlinearity was presented in |34j . The truncated equa- 
tions combine Eqs. ([t]) and Eqs. (|8]): 

dt^p = T(p,q,p_)v|/qVl/p e'^-* 

+ T(p,-q,p+)^qVl'p_e^^+* 

dt^^ = T(q,p,-p_)vl'p^p_e-^^-* 

-t- T(q,-p,p+)^pv|/p^e^^+* (9) 

dt-^p_ = T(p_,p,-q)*p¥qe-*^-* 

a^vl/p^ = T(p+,p,q)vI/pVl/qe-^A+t 

Strictly speaking, the chosen four modes {ipo, ipq, ip+ ^^^d 
■0-) are coupled to further modes and do not form a 
closed system. Indeed, even the linear problem closes 
only if all the satellites ±q -I- mp (m is a positive or 
negative integer) are included [1]. However, in consider- 
ing the linear instability it is traditional to truncate the 
system to the four modes only with a justification that 
the higher order satellites are less excited in the linear 
eigenvectors, which turns out to be a very good approx- 
imation if M ^ 1 and quite reasonable for M ^ 1 and 
larger [T]. In this paper we will test predictions of the 
4MT system, both linear and nonlinear, against DNS of 
the full system. 



4 



Nonlinearity parameter M 

In studying instability of a primary monochromatic 
wave, we will follow the convention that the wavenumber 
of this wave is p and denoted its amplitude by ^'p. The 
character of the instability is greatly influenced by the 
initial amplitude of the primary wave, ^'o = ^'p|t=o [I]- 
Following Gill [TI, we introduce the nonlinearity param- 
eter 



M 



5^ 



(10) 



M measures the relative strength of the linear and non- 
linear terms at the scale of the carrier wave. M ^ 1 
corresponds to the case where the /^-effect plays no role. 
For F — Q this case reduces to the Euler equation limit 
and the well-studied instability of the plane parallel sinu- 
soidal shear flow known as Kolmogorov flow [26 . Note 
that most papers on the modulational instability within 
the plasma context have dealt only with this limit (e.g. 
[20j EI])- Case M <C 1 corresponds to the weak non- 
linearity limit dominated by resonant wave triads. In 
this case the four constituent modes (carrier wave, mod- 
ulation and two satellites) can be split into two coupled 
triads which produce independent contributions to the 
instability |T]. The instability associated with a single 
triad is known as the decay instability |27]. The con- 
dition M ^ 1 deflnes the Rhines scale kr which marks 
a crossover from the hydrodynamic vortex to the wave 
behavior |35j. 

DECAY INSTABILITY OF A ROSSBY WAVE 

The decay instability is an instability of a carrier wave 
involving a pair of other modes (i.e. the primary wave 
decays into two secondary waves, see e.g. ^Z]). We 
shall derive this instability from the 3MT, Eqs. Q. In- 
troducing the vector notation ^ — (^fp, ^fq, ^'p_ ), a 
monochromatic carrier wave is given by 'S'o = (^'o,0, 0) 
where ^'g is a complex constant representing the ampli- 
tude of the initial carrier wave. This is an exact solution 
of Eqs. Q. We consider the stability of this solution 
to small perturbations involving the modes q and p_ 
by taking '3/ ~ 'U'o + e'3/i with the perturbation given 
by ^1 = (0, i/'q, ^/ip ). Linearisation yields the following 
equations at flrst order in e: 



i A_ t 



dtipq = T(q,p, -p_)*o'0p_e" 
ajp_ = T(p_,p,-q)^oV^qe^^-*. 
We now seek harmonic solutions: 



(11) 



V'p-(i) 



Ap_e 



This requires fip 



A_. Solving Eqs. (11) then 



reduces to flnding solutions of the linear system 

A, 



where 

A = 



= 



r(q,p,-p_)*o 



T(p_,p,-q)«'o i(-r!q + A_) 



(12) 



To obtain non-trivial solutions, we require detA — 0, 
which yields the dispersion relation: 

r!q(-r!q + A_) - T(q, p, -p_) T(p_, p, -q) l^ol" - 0. 

(13) 

This has two roots, ilq with corresponding eigenvectors: 



A, 



/ 1 

= T(p_^^p,^q) -^0 

\ i(Oq-A_) 



(14) 



Instability occurs when fiq has a non-zero imaginary 
part. For an exactly resonant triad, A_ = 0. For reso- 
nant triads, using Eq. ^ the roots of Eq. ([l3| are 



l*o| |P X q| 



(<Z2 + F)(p^ +F) 



(15) 

In this case, instability occurs if q < p < 

Before investigating the non-resonant instability fur- 
ther, it is convenient to perform some rescalings. The 
dimensionless carrier wave amplitude will be given by M 



defined in Eq. (|10|). We non-dimensionalise the other 



terms in Eq. (13) as follows: 



F 

p 
q 



pp, 



-n, 
p 

p^F, 



spq, 



where p = (px^Py) and q — {qx,Qy) are unit vec- 
tors pointing in the directions of p and q respectively. 



Eq. (13) can then be re-arranged to the following form: 
f^(-f^+A_)-M2r(sq,p,-p_)T(p_,p,-q) = 0. (16) 



where p_ 
roots are 



= p — sq and A_ — uj-p 



The two 



1 



A_ ± 



(A_)2 - 4M2T(sq,p, -p_) r(p_, p, -q) 

(17) 



To have an instability we require 

A- < 2M VT(sq, p, -p_ ) T(p_ , p, -q) 

which demonstrates that the instability concentrates on 
the resonant manifold, A_ = as M ^ 0. This is illus- 
trated in Fig. [l] The corresponding analysis for the triad 
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FIG. 1: The growth rate of the decay instabihty (the negative imaginary part of the roots of Eq.( 16 1) is plotted as a function of 
q for a fixed meridional carrier wave-vector, p = (1, 0), for various values of the nonlinearity parameter M. The set of unstable 
perturbations become concentrated on the resonant manifolds (blue lines) as the nonlinearity of the carrier wave is decreased. 



(p, — q, p+) produces identical surfaces reflected about 
the vertical axis reflecting the instability concentrating 
on the second resonant manifold, A+ = 0. As M — > 0, 
these two surfaces become disjoint from each other except 
near the origin q = 0. 



MODULATIONAL INSTABILITY: LINEAR 
ANALYSIS 

Let us now derive the modulational instability in the 
same way as we have done for the decay instability. The 
modulational instability is studied using the 4MT. We 
begin by linearising Eqs. (|9| about the pure carrier wave 
solution, = (^fQ, 0,0,0) where VPg is a complex con- 
stant representing the amplitude of the initial carrier 
wave. We consider the stability of this solution to small 



perturbations involving the 3 modes q, p_ and p+ by 
taking ^ = \I'o + e'^'i with the perturbation given by 
— (0, V'q, V'p4- ! V'p- )• Linearisation yields the follow- 
ing equations at first order in e: 

= r(q,p,-p_)M/oVe"'^-* 

+r(q,-p,p+)^7Ap^e^^+* (18) 

dt^p^ = r(p+,p,q)vi/o^qe-'^+* 

5jp_ = r(p_,p,-q)^o^qe^'^-*. 
We again seek harmonic solutions: 

^p,(<) = Ap^e-^V* 
V'pjt) = Ap„e-^^>'-*. 
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FIG. 2: Growth rate of the modulational instabihty (the negative imaginary part of the roots of Eq.(23l) as a function of q 
for a fixed meridional carrier wave-vector, p = (1, 0) and F — 0. The values of the nonlinearity M for the initial carrier wave 
are M = 10 (Euler limit), M = 1, M = 3/4, M = 1/2, A4 = 1/4 and M = 1/10 (weakly nonlinear limit). The set of unstable 
perturbations become concentrated on the resonant manifolds as the nonlinearity of the carrier wave is decreased. 
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This requires requires fip^ = riq + A+ and rip_ 



-T(q,-p,p+)T(p+,p,q) |*o| 



A_. Solving Eqs. (18 1 then reduces to finding solutions 
of the linear system 



A I A 



p+ 



= 



-T(q,p,-p_)r(p_,p,-q) |«'or(f^q + 
The corresponding eigenvectors are given by 



where 



A, 



1 

r(p+,p,q)1'o 
-i(flc|+A+) 
T(p_,p -q)^n 



(21) 



A 



T(p-,p,-q) ^-0 



i(f2q-A_) 



T{q,-p,p+)^o r(q,p,-p_)^'o 
iin^ + A+) 

-i(-r2q + A_) 

(19) This derivation holds for any system with a quadratic 
Setting det ^ = yields a cubic dispersion relation: nonlinearity. Using Eq. ^ and performing some alge- 



nq(nq + A + )(-r!q + A_) 



bra we recover the usual form of the dispersion relation 
(20) specific to the CHM equation P (see also pi IT51 [2CT1 1^ ) : 



(g2 + F)n + f3q., + \^o\ |P X qr (p' - q^) 



{pl + F){n + Lo) + pp+^ {pi + F){n - uj) + I3p^ 



= 



(22) 



r 



This can be solved numerically, and sometimes analyti- 
cally, for a given set of parameters to determine fi. For 
the purposes of easy comparison of diff'erent values of M, 
we nondimensionalise as before. The result is 

(s^ + F)n + sq^ + M'^s'^il - s^) \p X q|^ [r+ - T"] = 0, 

(23) 

where 



|p±sqr-l 



(|p ± sq| + F){-^ ±n)+p,± sq.. 



(24) 



The roots of this equation are controlled by five param- 
eters, M, F, s, 6p and 9q where 9p and 6q are the angles 
between the x-axis and the carrier wave-vector and per- 
turbation wave-vector respectively. The structure of the 
instability is strongly dependent on the value of M. This 
is shown in Fig. |2] We see that the modulational insta- 
bility is, in some sense, the nonlinear sum of the decay 
instabilities for the two triads, and we will clarify this 
issue in the next section. 



COMPARISON OF THE 3MT AND THE 4MT 
MODELS WITH DNS OF THE FULL CHM 
SYSTEM 

There is sometimes confusion in the literature, perhaps 
partially semantic, on whether the modulational instabil- 
ity of the Rossby and drift waves is governed by 3-wave 
or 4- wave interactions. Here we will clarify this issue. It 
was shown by Gill that as M ^ 0, the modulational 



instability obtained within the 4MT model localises on 
the two resonant manifolds for the two triads A+ = and 
A_ = 0. Since these two curves are mostly disjoint from 
each other (except for the origin), in the weakly nonlinear 
limit, the modulational instability is just a simple sum 
of the two decay instabilities. Namely, the two unstable 
eigenvectors of the instability of the 4MT will coincide 
with the eigenvectors of the two respective branches of 
the decay instability (i.e. the fourth mode in such 4MT 
eigenvectors will have zero amplitude). In particular, the 
maximum growth rates of the 3MT and 4MT instabilities 
become identical. For larger values of M, the growth rate 
of the modulational instability is typically larger than 
that of the corresponding decay instability. 

However, for the typical setup where the primary wave 
is purely meridional and the modulation is purely zonal, 
the wavevector q is equally close to both branches of 
the three-wave resonant manifold. This is because these 
resonant manifolds cross zero of the q-space in the di- 
rection of the gj,-axis, i.e. in the zonal direction. Thus, 
the above speculations about the equivalence of the 3MT 
and the 4MT for weak waves may not apply to such a 
setup. Therefore, let us consider the weakly nonlinear 
case (M = 0.1) and examine predictions of the 3MT and 
the 4MT models and compare them to DNS of the full 
CHM system in the following two cases: 

(A) the primary wave is purely meridional, p — (10, 0), 
and the modulation is purely zonal, q = (0, 1); and 

(B) the primary wave is purely meridional, p = (10, 0), 
and the modulation is off- zonal. We take q = (9, 6). 
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FIG. 3: Amplitude of the zonal mode with wavenumber q for 
M = 0.1, p = (10,0) obtained from DNS and from solutions 
of 3MT and 4MT models. Case (i): purely zonal modulations, 
q=(0,l). 
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FIG. 4: Comparison of modulational and decay instabilities 
for early times for the case (i) -purely zonal modulations (like 
in Fig.[3|. 



This is close to the maximum of the most unstable 
mode on the resonant curve We cannot select the 
exact value of the maximum because the discrete 
wave numbers in the periodic box do not typically 
lie exactly on the resonant manifolds. This is a sub- 
tlety which can have strong implications for numer- 
ical simulations of very weakly nonlinear regimes 
|36j which we have been careful to avoid here. 

Let us first consider case (A) when the modulation is 
purely zonal. Fig. |3] compares |\['q| obtained from the 
solution of Eq. ([3| with that obtained from solutions of 
3MT, Eqs. (0, and 4MT, Eqs. (|9|. The initial condition 



□ IVI = 0.1 4-mode truncation 
IVI = 0.1 3-mode truncation 
IVI = 0.1 (DNS) 




FIG. 5: Same as in Fig. [3] but now for the case (ii): off-zonal 
modulations, q = (9, 6). 



was constructed from Eq. ( 14 1 , the unstable eigenvector 



for the decay instability. We see that the growth rate pre- 
dicted for the decay instability is not observed. The PDE 
instead seems to follow the growth rate for the modula- 
tional instability. From the zoomed-in plot of the early 
time evolution shown in Fig. |4] we see that the full dy- 
namics very quickly generates the mode p-)_ which is ab- 
sent from Eqs. ^ . The full system then quickly deviates 
from the solution of the 3MT in a time of the order of 
the inverse of the instability growth rate. However, the 
set of 4 modes takes much longer to generate any further 
modes. Thus in this setup, the 4MT Eqs. (|9| provide a 
much better description of the full dynamics for times up 
to 10 instability times. 

Let us now consider the case (B) when the modula- 
tion is off-zonal. Fig. |5] compares |^'q| obtained from the 
solution of Eq. (|3| with that obtained from solutions of 
3MT, Eqs. (0, and 4MT, Eqs. for an initial condi- 
tion being the unstable eigenvector for the decay insta- 
bility, Eq. ([U]). As expected, now the 3MT and the 4MT 
give practically identical results, and both of these mod- 
els agree well with DNS up to the time equal to seven 
inverse growthrates. They predict well the maximum of 
the zonal jet amplitude, although the subsequent stage 
of decrease is not described as well as in case (A) by the 
4MT model. 

From these results we conclude that the 3-wave inter- 
action is indeed the basic nonlinear process when M <C 1 
provided the triad is not degenerate, in the sense that it 
does not contain quasi-resonant modes which are equidis- 
tant from two different resonant manifolds as happens 
when the vector q is zonal. In these cases, the 3MT 
system is just as good as the 4MT and it describes well 
the full CHM system for over several characteristic times 
(i.e. the inverse instability growthrates). On the other 
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hand, the most relevant configuration with q zonal is, in 
fact, degenerate. In this case, however, the 4MT model 
works well over many characteristic times whereas the 
3MT fails almost immediately. Thus, to have a wider 
range of applicability, we will study the 4MT model and 
abandon the 3MT model in the remainder of the present 
paper. 

Next, let us study the modulational instability arising 
from the 4MT in greater detail. 



INSTABILITY FOR PURELY MERIDIONAL 
CARRIER WAVE AND PURELY ZONAL 
MODULATION 



The case of a purely zonal carrier wave (p = (1,0)) and 
purely meridional perturbation(q = (0, 1)) is of physical 
interest and produces considerable simplification. The 
dispersion relation then reduces to solving 



which has roots 



(1 + F)(s2 + 1 + F) 



2AP{1 - s2)(l + F)2(s2 + F + 1) - (s2 + F) 



n = o. 



n = ±i 



(1 + F)(s2 + 1 + F) 



2M2(1 - s2)(l + F)2(s2 + F + 1) - (s2 + F) 

s2 + F 



(25) 

(26) 
(27) 



The question of whether the perturbation is unstable reduces to the question of when the quantity under the square 
root is positive. In this expression, recall that s is the ratio, q/p, of the modulus of the modulation wave-vector to 
the modulus of the primary wave-vector. Letting s'^ = y, one obtains a quadratic for the quantity under the square 
root which is positive in the range s G (— Smaxi Smax) where 



1 -I- 2M2F(1 + F)2 
2M2(1-KF)2 




(2A/2(1 + Ff - F)(2M2(1 + FY 
(H-2M2F(H-F)2)2 



(28) 



Case of Infinite Deformation Radius 

When F = the analysis becomes particularly simple. 
There is always a range of unstable long wavelength per- 
turbations, given by (0, Smax), for any value of M . Smax 
is given by 



/-I + Vl + 16M4 
Within this range the growth rate is 



n = 



(s2 + 



v/2M2(l - s4) 



The growth rate has a single maximum at sq 
where ?/o is the positive root of 

2/3 + 32/2 + (1 + -L)y_ 1=0. 



(29) 
(30) 



(31) 



One can show that so^vv2— lasM— >oo and sq = 
M + 0{M^) as M ^ 0. One would be interested to know 
when the maximally unstable meridional perturbation is 
a local maximum with respect to nearby non-meridional 



perturbations. To ascertain this, one should look at the 
sign of the determinant 



AM{qx,qy) 



dq^dgy d^tjy 



(32) 



evaluated at {qx,qy) — (0,so). This can be done semi- 
analytically using Mathematica and is plotted in the inset 
of Fig. [e] We find that Am > with ^ < (the 
criterion for a local maximum) for M > M^. Am < 
with < (the criterion for a saddle) for M > 

Mc- The critical value of M is found numerically to be 
Mc w 0.534734. Numerical explorations show that the 
local maximum found for M > Mc, is actually global. 
For M > Mc, therefore, the fastest growing perturbation 
is indeed zonal. As M decreases below Mc the most 
unstable pertubation moves to a point with a finite value 
of qx- The maximally unstable perturbation for M < 
0.53 tends to a point on the resonant manifold making 
an angle of 57r/6 with the x-axis. The dependence of this 
angle on M is shown in Fig. [6j A clear transition from an 
axial maximum to an off-axis maximum is clearly visible. 
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71/3 - 
71/6 - 

I ' ' ' ' 

0.2 0.4 0.6 0.8 1 

M 

FIG. 6: Angle, 9, between the q wave- vector of the maxi- 
mally unstable perturbation and the a;-axis as a function of 
M. Inset plots Am and flxx as a function of M illustrating 
the transition of the maximum growth rate for on-axis pertur- 
bations from a local maximum to a saddle point at M « 0.53. 



F=OA — B- 
F=0.5 —e- 




0.2 0.4 0.6 0.8 1 1.2 



s 

FIG. 7: Instability growthrate for purely meridional pertur- 
bations with M = 0.5 > •\/2/27 for different values of the 
deformation radius. 



Effects of Finite Deformation Radius 

We now consider the dependence of MI on the deforma- 
tion or Larmor radius, noting that a finite deformation 
radius is obtained in the QG system under a reduced 
gravity approximation. When F is finite, there are 2 
regimes, depending on the value of AI . For an interval 
of instability to exist, we require s'^-^^^ > 0. This requires 
that 



p{F) = 2M2(1 + Ff -F >0. 



(33) 



The discriminant of the corresponding cubic, p{F) = 0, 
is — 4(— 2M^ -I- 27M'*) . Since we are only interested in 
> we can identify two regimes. 



p{F) — has one real root, 
> when F > 



(281 



Regime 1: M > 
Referring to Eq 

Fi (which is negative) and p{F) 
Fi. Then for any positive value of F there exists 
a finite range of s, s G (0,Sniax)i for which the 
perturbation is unstable. Smax is given by Eq. (28 1. 
In this regime, finite deformation radius tends to 
reduce the growth rate of the instability but cannot 
suppress it. See Fig. [7] 

Regime 2: M < 
Referring to Eq. ( |28[ ), p{F) = has three real roots, 
Fi, F2 and F3. Fi is negative and F2 and F3 are 
positive. p{F) < in the range (^2,^3). In this 
regime, there are critical values of F, Fi and F2 
such that the range s S (0, Smax) of unstable pertu- 
bations only exists if F < Fi or F > F2 . Fi and F2 



0.06 



0.05 



0.04 - 



■ 0.03 



0.02 



0.01 




FIG. 8: Instability growthrate for purely meridional per- 
turbations with M = 0.25 < •y/2/27 for different values of 
the deformation radius. For this value of M, F2 ~ 0.23 and 
F3 ~ 1.00. Note that the instability is completely suppressed 
for intermediate values of F and then emerges again as F 
increases. 



are obtained by finding the positive roots of Eq. 33 
and Smax is again given by Eq. (28). In this regime. 



there is a range of intermediate deformation radii 
which completely suppresses the instability. See 
Fig-i 




FIG. 9: Same as in Fig. |2] but now for a finite deformation radius, F = 2. 
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FIG. 10: Growth rate of the modulational instability (found from Eq.(23l) as a function of p for a fixed zonal modulation 
wave-vector, q = (0, 1) and F = 0. The values of M for the initial carrier wave are A'l = 10 (Euler limit), M = 1, A4" = 3/4, 
M = 1/2, M = 1/4 and M = 1/10 (weakly nonhnear hmit). 




FIG. 11: Same as in Fig. |10| but now for a finite deformation radius, F = 2. 
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ROLE OF THE CARRIER WAVE AMPLITUDE. 

We have already mentioned the role of the nonlinear- 
ity parameter M, for the most unstable modulation (i.e. 
zonal for big M and inclined for small M), as well as 
for the different regimes in the finite F case. Fig. [2] 
show plots of the instability growthrate as a function of 
1 = {Qx,qy) for several different values of M for fixed 
meridional p and F = 0. In particular we see how the 
maximum growthrate flips from zonal to off-zonal q when 
M is reduced (below see also about the collapse of the 
unstable region to the resonant curve). Fig.[9]shows plots 
of the instability growthrate similar to the ones in Fig. [2] 
but now for finite F. Qualitatively the finite F plots are 
similar to those obtained for F = 0. 

Another natural way to visualise the structure of the 
set of unstable pertubations is to fix the wavevector of 
the perturbation mode, q, and plot the instability growth 
rate as a function of the primary wavevector, p. Fig. [TO] 
does this, plotting the instability growth rate as a func- 
tion of p = {px,Py) for several different values of M with 
fixed zonal q and F — 0. We see that nonlinearity re- 
duction "eats into" the instability cone, i.e. makes some 
wavenumbers inside the cone stable. At the same time, 
the nonlinearity makes some wavenumbers outside the 
cone unstable. It is important to keep in mind that, 
even for large M, the maximum growthrate occurs out- 
side of the cone, for the primary wave orientations closer 
to zonal than to the the meridional direction, see Fig. [TO] 
for M = 10. This fact is easy to overlook if one considers 
only the limit M — > oo (as it is common in the plasma lit- 
erature) because, in this limit, the growth rate maximum 
is for the meridional primary waves. On the other hand, 
the choice of the primary wave direction is often dictated 
not by the maximum growthrate of the modulational in- 
stability, but by the structure of the primary instability 
creating the Rossy and drift waves (ITG instability in 
plasmas and the baroclinic instability in GFD). 

Finally, Fig. [TT] shows plots of the instability growth 
rate similar to the ones in Fig. [TO] but now for finite F. 
We again see a qualitatively similar picture to the F = 
case. Note, however, that the global maximum growth 
rate for large M is now obtained for purely meridional 
primary waves. 

Let us now specially consider the limits M ^ 1 and 
M < 1. 

Limit M > 1. 

The limit of large nonlinearity M ^ 1 is a particularly 
simple and well studied one [II [HI |20l llll 126]. As 

we mentioned before, the /3-effect becomes unimportant 
and, for F — 0, this case reduces to instability of Kol- 
mogorov flow in Euler equations (i.e. sinusoidal plane- 
parallel shear). In this case the most unstable modula- 



tion is perpendicular to the carrier wave. The instability 
criterion reduces to [1] 

cos2</>< (^1 + ^) /4, 

where (j) is the angle between p and q. For the scale sepa- 
rated case, q <^p, this condition describes an "instabihty 
cone" [2l|20l[2T] 

101 <^/6. 

Finite deformation radius modifles this cone to a larger 
instability area [2 

F + pI~3pI>0. 

We repeat that one has to use the results obtained in 
the limit M ^ oo with great caution, because even for 
large M's the most unstable primary wave is not pre- 
dicted correctly in this limit. 

Limit Af < 1. 

In the limit of weak nonlinearity, Af <C 1, the dy- 
namics is completely wave dominated [1]. The nonlinear 
terms allow waves to interact weakly and exchange en- 
ergy. Since the nonlinearity is quadratic, wave interac- 
tions are triadic (3-wave resonances are allowed by the 
dispersion relation, Eq. Q). Any triad of waves having 
wave- vectors ki, k2 and ks interact only if they satisfy 
the resonance conditions: 

k3 = ki + k2 (34) 

w(k3) = w(ki)+a;(k2). (35) 

From Eq.(|4|, this latter relation gives an implicit equa- 
tion for the resonant manifold of a given ka — (fca^, k^y): 

kix ^ kix 
kfx + kly + F {kax - hx)^ + {ksy ~ kly)^ + F 

-T^^tI S = 0- (36) 

klx + kly + F 

Because the system is anisotropic, the shape of resonant 
manifold depends on the direction of ka as shown in 
Fig. [121 

These resonant manifolds are relevant even for flnite 
nonlinearity since the support of the instability concen- 
trates close to the resonant curves as M is decreased as 
shown in Fig. [2] Even for M = 1 there is a strong connec- 
tion between the resonant curves and the shape of the set 
of modulationally unstable perturbations. In fact. Fig. [2] 
shows two resonant curves corresponding to two resonant 
triads, 

ki = p, k2 = q and ks = p+ = p -I- q 
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This fact can be explained by considering the resonant 



curve for q <^ p where it behaves as Qx 



-2{pxPy/p^)ql 



For instabihty, this curve has to pass as close as possi- 
ble to the vertical (zonal) axis (where we have chosen 
our q). Thus, we need to minimize the above coefficient 
{PxPy/p'^) (e.g. with respect to py for fixed Px) which 
immediately gives Py — ipx/VS. 

For small M the maximally unstable modulation q is 
off-zonal, which may be important for determining the fi- 
nal statistical state of the nonlinear evolution. As we will 
see later, this state appears to have a predominantly off- 
zonal component even if the initial modulation is chosen 
to be zonal. 



NONLINEAR EVOLUTION 



FIG. 12: Shape of the resonant manifold determined by 
Eq.(36l with ka — (cos(&), sin(^)) for several values of 9 for 
the case F — 0. 



and 



ki 



and 



Out of four wavenumbers in our truncated system, 
PiQiP- E^J^d p_|_, three are resonant (or nearly reso- 
nant) and the remaining one is non-resonant (p_ or p+). 
As we mentioned before, this picture is correct in non- 
degenerate situations, when q is not zonal. Then for 
M the amplitude of this non-rcsonant mode in the 
instability eigenvector tends to zero, so effectively there 
are only three active modes, and one can use the re- 
sults obtained above for the 3MT model. In particular. 



Eq. (15) gives the instability growthrate: 



IV'ollq X p\^J {p^ ~ q^){pl - 
^ipl + F){q^ + F) 



(37) 



This expression was previously obtained in the case F — 
in [11 based on the 4MT model. One can see that 
instability of the primary wave occurs if its wavenumber 
length lies in between of the wavenumber lengths of the 
waves it decays into, q < p < p-^-. This condition has 
a nice dual-cascade interpretation: to decay the wave 
must be able to transfer its energy to a large scale and 
its enstrophy to a smaller scale. For F = 0, the typical 
instability growthrate is 7 ~ UqP where Uq = pipo is the 
velocity amplitude of the carrier wave [T]. In the large 
F case, the instability is slowed by the factor F/p'^ (but 
not arrested). 

Another interesting feature of instability for Af ^ 1 is 
seen in Fig [To| where we can see that (for fixed zonal q): 
the unstable region becomes narrow and collapses onto 
the sides of the "cone", i.e. onto the lines Py = ±Px/V^- 



From now on we study only systems with infinite defor- 
mation radius, F = 0. To test the linear predictions, and 
to study the nonlinear evolution, we have performed DNS 
of the CHM system, Eq. ([l}, using a standard pseudo- 
spectral method with resolution up to 1024^ and hyper- 
viscosity parameters i^n = 4.5e~^°. We solve, in tan- 
dem, the 4MT system, and compare it with DNS. 
Although the 4MT was used as the departure point for 
the linear stability analysis, it is a fully nonlinear set of 
equations in its own right. In addition to checking the 
linear instability predictions against DNS, we will also 
explore the extent to which the nonlinear dynamics of 
the 4MT captures the behaviour of the full PDE. In all 
cases, we choose the initial condition to be along the un- 
stable eigenvector of the 4MT. 



Case of Meridional Carrier Wave and Zonal 
Modulation. 

Let us first of all consider the geometry which we 
dealt with most: purely meridional carrier wave and 
purely zonal modulation. We choose p = (10, 0) and 
q = (0, 1). A series of frames of the vorticity field for the 
cases of strong {M — 10), medium [M = 1) and weak 
(M = 0.1) nonlinearities obtained by DNS are shown in 
Figs. (13 1, (14 1 and (151 respectively. The evolution of 



the mean zonal velocity u{y), averaged over x, obtained 
from DNS for the same set of nonlinearities is shown in 
Figs. (16 1, (17 1 and (18 1 respectively for times close to 



the formation of the jet. Finally, evolution of the ampli- 
tude, iV'ql, of the zonal mode for the same runs is shown 



in Figs (19), (20) and (21) respectively. For compari 



son, we also put the corresponding values of IV'ql obtained 
from the 4MT. 

Immediately, one can see that the initial stage of evo- 
lution agrees very well withe predictions of the linear the- 
ory obtained from the 4MT. Moreover, the 4MT works 
rather well beyond the linear stage, particularly in the 




FIG. 13: Vorticity snapshots showing the growth, saturation and transition to turbulence of a zonal perturbation of a meridional 
carrier wave having M = 10. 




FIG. 14: Vorticity snapshots showing the growth, saturation and transition to turbulence of a zonal perturbation of a meridional 
carrier wave having M = 1. 
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FIG. 15: Vorticity snapshots showing the growth, saturation and transition to turbulence of a zonal perturbation of a meridional 
carrier wave having Af = 0.1. 
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FIG. 16: Mean zonal velocity for M = 10 



FIG. 18: Mean zonal velocity for M = 0.1 
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FIG. 17: Mean zonal velocity for M = 1 



M = 0.1 case, where the initial growth reverses in agree- 
ment with the (periodic) behavior of the four-mode sys- 
tem. For M = 1, the system's growth does not reverse, 
but rather experiences a saturation at the level where 
the four-mode system reaches maximum and reverses. 
The most surprising behavior is seen for M = 10 where 
the linear exponential growth continues well beyond the 
point of reversal of the four-wave system, even though 
the system is clearly nonlinear at these times and follows 
a self-similar evolution, see below. 

A common feature of the nonlinear saturation stage 
of the jet growth is self-focusing of the zonal jets which 
become very narrow with respect to the initial modula- 
tion wavelength. This self-focusing cannot be described 
by the 4MT because such anharmonic jet shapes involve 
strong contributions from higher harmonics p ± nq. For 



M = 10 4-mode truncation 
M = 10 DNS 




FIG. 19: Growth of the zonal mode q obtained by DNS and 
by solving 4MT system for p = (10,0) and q = (0,1) and 
for M — 10. Time has been scaled by r (the inverse of the 
instability growth rate). 



large M and q <^ k, such jet "pinching" was predicted 
theoretically in [5] where self-similar solutions were ob- 
tained describing a collapse of the jet width. These 
strong narrow zonal jets are expected to produce trans- 
port barriers in the transverse (y) direction, which is im- 
portant in both fusion plasma and the geophysical con- 
texts. 



Figure ( 22 1 shows the zonal velocity u re-scaled with 



self-similar variables as u{y,t) — a(t) f{b{t)y) in the run 
with M = 10. The self-similar stage occurs in the 



time interval corresponding to the overshoot in Fig (19 1, 
i.e. after the 4MT has reached its maximum but be- 
fore DNS saturated at a plateau. Empirically, we obtain 

At least during the decade 



a{t) = uq e'^* and b{t) 



,1.85* 
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M = 1 4-mode truncation 
M = 1 DNS 
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FIG. 20; Same as in Fig. [19] but for M = 1 
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FIG. 21: Same as in Fig. [19] but for M = 0.1. DNS results 
are presented from calculations at two different resolutions to 
illustrate that the oscillatory dynamics of the zonal mode are 
not influenced by small scale dissipation. 



of amplitude growth observed before pinching occurs, one 
can see a good evidence of self-similar behaviour and, re- 
markably, the nonlinear growth at the self-similar stage 
follows the same exponential law with growthrate 7 as on 
the linear instability stage. Note that the self-similar so- 
lutions were obtained in [2] based on the scale separated 
description and, therefore, the self-similar pinching must 
stop when the scale separation property breaks down due 
to the jet narrowing (at which point a roll-up into vortices 
occurs, see below). In the smaller M runs, the overshoot 
is absent and the amplitude of the zonal mode decreases 
after reaching a maximum in correspondence with the 
solution of the 4MT. The self- focusing is thereby much 
reduced and the self-similar stage is not clearly observed. 
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FIG. 22: Zonal velocity u(y) re-scaled with self-similar vari- 
ables: u{y,t) = a{t) f{b{t)y) with a{t) = uoe^* and b{t) = 
e^-^^^ for M = 10. 



One can also see some qualitative behavior differences 
for different degrees of nonlinearity M. First of all, we see 
that the east-west asymmetry is larger for weaker waves, 
which is seen as asymmetry of the top and bottom halves 
of the vorticity distributions in Figs. 13] 14 and 15 This 
is natural considering that for large nonlinearities the 
beta-effect, which is the cause of the east- west asymme- 
try, is less important. Further, we see that for large M 
the nonlinear evolution is vortex dominated and that the 
vorticity of the initial carrier wave rolls into vortices orga- 
nized into Karman-like vortex streets. This corresponds 
to the moment when the jet velocity reaches a plateau 
in Figs (19) and (20 1. On the other hand, in the weak 



wave case M <C 1 one cannot see vortex roU-ups and the 
dynamics remains wave dominated. 

For large nonlinearities, at the final stages the vortex 
streets break up due to a vortex pairing instability, which 
is followed by a transition to turbulence. Such turbulence 
is anisotropic with a pronounced zonal jet component. 
On the late frames Fig. [13] we can see a well formed po- 
tential vorticity staircase structure as described in |37j . 

For small nonlinearities, M <C 1, the system's non- 
linear evolution starts with self-focusing of the primary 
wave, but this is followed by a quasi-oscillatory behav- 
ior where the system returns close to the initial state. 
This is very well reproduced by the four-wave trunca- 
tion. The same effect was also noted for the Generalised 
Hasegawa-Mima model in [SI] and in the atmospheric 
dynamics context in |38| . However, the periodic behav- 
ior is not sustained and a transition to an anisotropic 
turbulent state occurs. It is interesting that the domi- 
nant jet structures observed in such a turbulent state are 
off-zonal. This effect may be connected to the off-zonal 
"striations" reported for the ocean observations in [7]. 
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We currently regard this connection with some caution 
since we have not performed any time averaging whereas 
the ocean striations are sufficiently weak that they only 
become evident in the averaged data. We hope to inves- 
tigate this further in future work. 

For Af > 1, the vortex streets represent the 2D fine 
structure of the saturated zonal jets (i.e. at the plateau 
part of Figs (19 1 and (20l). Such vortex street con- 
figurations are more stable than the plane parallel {x- 
independent) flows with the same zonal velocity proflles. 
This can be understood heuristically (see, for example 
[55] , chap. 3) by considering the vortices to impart some 
eddy viscosity to the mean zonal velocity profile. Recall 
that stability of the latter is determined by the Rayleigh- 
Kuo necessary instability condition [SH], 



dyyU{y) -P>0. 



(38) 



Figs (241, (251 and (261 plot the profiles of dyyu{y) ~ f3 



at different moments in time corresponding to runs with 
M = 10, M = 1 and M = 0.1 respectively. One can 



see that in Figs (24) and (25) these profiles cross the x- 
axis (especially far in the M = 10 case) which implies 
that the zonal flows get stronger than the limiting val- 
ues implied by the Rayleigh-Kuo condition. We interpret 
this as a result of a competition between the instability 
and the jet pinching process. For large M the pinching 
is self-accelerating (self-similar) and it manages to sig- 
niflcantly compress/ amplify the unstable jet in the finite 
time needed for the instability to develop (i.e. the inverse 
growthrate). On the other hand, in the case M = 0.1 the 
jet strength reaches a maximum and then decreases re- 
maining in the stable range according to the criterion 
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These results allow us to draw conclusions about the 
critical value of nonlinearity, M = , which separates 
the two qualitatively different types of behavior: vor- 
tex roll-up and saturation vs the oscillatory dynamics, 
see Figure (23). If the jet strength maximum, as pre- 



dicted by the 4MT, exceeds the values of Rayleigh-Kuo 



necessary instability condition ( 38 1 then the vortex roll 



up occurs and the jet enters into a saturated, relatively 
long-lived plateau stage. At this moment, the system's 
behavior starts to depart from the 4MT model. On the 
other hand, if the jet strength maximum, as predicted 
by the 4MT, remains below the Rayleigh-Kuo threshold 
then the system's growth reverses and it follows the 4MT 
dynamics for longer time. 

This simple picture leads to a qualitative physical es- 
timate for and for the saturated velocity of the jet. 



Let us start with the latter, see Fig. (23) 



Since the x-periodicity is preserved, the step of the 
vortex street is equal to the wavelength of the original 
carrier wave. The vortices in the stable vortex street are 
approximately round and the y spacing between the vor- 
tices is approximately the same as the x-spacing. Thus, 
the saturation width of the pinched jet is of order of the 
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FIG. 23: Growth of zonal perturbations due to modulational 
instability of a meridional carrier wave having p = (10,0) 
for several different values of M. The amplitude of the zonal 
mode has been scaled by ^to and time has been scaled by r 
(the inverse of the instability growth rate). 




FIG. 24: The Rayleigh-Kuo profiles for M = 10 



wavelength of the initial carrier wave. Lagrangian con- 
servation of the potential vorticity determines the final 
saturated amplitude of the jet. Indeed, the same vortic- 
ity as in the initial carrier wave rolls into the vortices (the 
(3y part of the potential vorticity does not play much role 
here since the fluid parcels remain at about the same y's) 
and the rest of the vorticity is shed in between of the vor- 
tex streets, shredded by shearing and dissipated. Thus, 
the jet saturation velocity is of the order of the velocity 
amplitude of the initial carrier wave. 
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FIG. 25: The Rayleigh-Kuo profiles for M = 1 
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FIG. 26: The Rayleigh-Kuo profiles for M = 0.1 

This estimate is well confirmed by our numerical results 
for M = 1 and M = 10. Indeed, taking values of Umax 



from Figs (16 1 and (17l) we get 

- 

Umax ~ O 9 ■ 



(40) 



Now, estimating dyyu{y) as p Umax and using (40) we 



can rewrite the instability condition ( 38 1 in a very simple 
form as 



M > Ah ' 

Our numerics show that A'h 



1/3. 



0.25 - 0.35, see Fig- 
ure (23 1. Note that the boundary is not sharp. For 



M — 0.25 the dynamics is definitely wave-dominated, 
however some elongated fuzzy vortices are still apparent 
whereas for M = 0.35 streets of round vortices are clearly 
formed with some wave-like oscillations still present. 



Case of Meridional Carrier Wave and OfF-Zonal 
Modulation. 



Above we considered the case when the carrier wave 
is purely meridional and the modulation is purely zonal. 
This geometry is important considering that both the 
baroclinic instability in GFD and the drift-wave insta- 
bilities in plasmas typically have most unstable modes 
being in the meridional direction. These modes can be 
considered as an initial condition for the secondary mod- 
ulational instabilities as it is done in the present paper. 
At the same time, we have established above that the 
most unstable modulations for M > 0.53 are zonal. 

On the other hand, for low M the most unstable 
modulations are off-zonal. This, in our opinion, is the 
reason why the final statistical state in the system in 
the M — 0.1 simulation showed presence of off-zonal 
anisotropic flows even though the initial modulation was 
purely zonal. Moreover, it is quite likely that in such 
weakly-nonlinear cases the system will pick the modula- 
tion which is off-zonal already at the initial moments. 

Thus, here we will consider a case with M = 0.1 where 
we start with purely meridional carrier wave, p = (10, 0) 
and with the modulation wavevector corresponding to 
the fastest growing mode in this case, namely q — (9, 6). 
Corresponding numerical results for this case are shown 
in Fig. (27) (vorticity snapshots) and Fig. ^ (evolution 
of the g-modc amplitude \^pq\ and respective results ob- 
tained from simulating the 4MT and 3MT models). 



First of all, as in all previous cases, we see good agree- 
ment of the initial evolution with predictions for the lin- 
ear instability obtained based on the 4MT and the 3MT 
models. Moreover, we see that the 4MT and the 3MT 
in this case qualitatively describe the nonlinear behav- 
ior too. Namely, like in the four-mode system, we see 
oscillatory behavior, even though the oscillations appear 
to be irregular. However, these irregular oscillations are 
clearly non-turbulent, as one can see from the vorticity 



frames in Fig. ( 27 1 which shows quite a regular pattern 



even at t ~ 100 (in the units of the inverse instability 
growthrate), - by which time the respective M = 0.1 sys- 
tem with zonal q is completely turbulent, see Fig. [15] 
Another way to see that the dynamics are regular in this 



case is to look at the 2D A;-spectra shown in Fig. (28). 
At i = 0, the only excited modes are the carrier wave 
p, modulation q and two satellites p ± q: these modes 
are marked by bold symbols in Fig ( 28 1 . At t = 60 one 



can see a regular "crystalline" structure corresponding 
to a discrete set of nodes np + mv (with integer values 
of m and n) with energy within 1% of the initial car- 
rier wave energy. Transition to turbulence does eventu- 
ally occur after a very long time, and the turbulent state 
does exhibit off-zonal striations similar to the respective 
M = 0.1 system with zonal initial modulations q. 




FIG. 27: Vorticity snapshots showing the growth, saturation and transition to turbulence of an off-zonal perturbation of a 
meridional carrier wave having M = 0.1. 



24 



40 
30 
20 
10 
' 
-10 
-20 
-30 
-40 



"T 1 1 1 1 1 r 



Xq Xp^ 
Xp 



_l I I I I I L_ 



5 10 15 20 25 30 35 40 



0.01 



0.001 




M = 10 stable DNS 
0.015exp°^'; 

10 15 
t/x 



20 



FIG. 29: Growth of the zonal mode q obtained by DNS for 
p = (8, 6) and q = (0, 1) and for M = 10. 



FIG. 28: Excited modes at f = 60 with energy within 1% of 
the energy of the initial perturbation in run with for q — (9, 6), 
M = 0.1. 



STABLE CASE 

Above we considered in detail various situations where 
the hnear theory based on the 4MT model predicts in- 
stability. We have also investigated the linearly stable 
case. 

For small M the zonal mode in the modulationally 
stable case behaves as expected, following the 4MT the- 
ory without growth of the mode. In this case, deviations 
from the 4MT are tiny, hence we omitted the correspond- 
ing graph. For M ^ 1, the situation is more interesting. 
Fig. 29 shows the evolution of the zonal mode for the 
i, 6) and q 



run with p = (8, 6) and q — (0, 1) and for M = 10 which 
corresponds to a linearly stable configuration within the 
4MT model. We see agreement with the 4MT stabil- 
ity prediction at early times, i.e. the zonal mode is not 
for t < 1. 



growing in Fig. 
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However, after about one 
timescale the zonal mode quickly breaks into growth, in- 
creasing (more or less exponentially) by two orders of 
magnitude Hence, the 4MT instability criterion must be 
used with caution if M 1. Further, for M ^ 1 stable 
case Manin and Nazarenko [2] predicted zonal velocity 
profile steepening for and this is evident in Fig. |30| where 
the initial sinusoidal profile develops into a triangular 
Burger's shock-type profile. 



SUMMARY AND CONCLUSIONS 

In this paper we dealt with the theory and numer- 
ical simulations of the modulational instability of the 
Rossby/drift waves described by the CHM model. We 
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FIG. 30: Mean zonal velocity profile for the stable configu- 
ration p = (8, 6), q = (0, 1) and M = 10. 



have revisited the linear theory of Gill [1 using the 4- 
mode truncation and emphasised the role of the carrier 
wave amplitude/nonlinearity, the role of the deformation 
radius and the role of resonant wave interactions in the 
case of weakly nonlinear carrier wave. We found a change 
of the most unstable modulation from zonal to off- zonal 
when the carrier wave nonlinearity parameter M falls 
below a critical value, M > 0.53. This latter effect may 
be important for understanding the recent ocean obser- 
vation of off- zonal jet striations [7 . It is also a likely 
mechanism for generation of off- zonal random jets in our 
numerical simulations at the late development stages of 
the modulational instability for M = 0.1 case. 

We established how the modulational instability re- 
lates to the decay instability obtained within the 3-mode 
truncation in order to clarify the question whether the 
dominant nonlinear mechanism of the modulational in- 
stability is three-wave or four-wave. The 3MT works very 
well for low nonlinearities M when the carrier wave and 
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the modulation belong to the same resonant triad which 
is non-degenerate, i.e. when it does not include wavevec- 
tors too close to fc = point where the two branches 
of the resonant curve intersect. This excludes the most 
popular choice of purely zonal modulations, for which 
3MT appears to be a bad model. On the other hand, 
4MT is more general and it works very well for small 
nonlinearities M including the case of the purely zonal 
modulations. Moreover, 4MT also works well for the ini- 
tial evolution in the strongly nonlinear cases, M > 1, 
including the linear growth stage and the prediction of 
the critical nonlinearity Af* for which transition from the 
saturated to the oscillatory nonlinear regimes is observed. 

We checked the solutions of the truncated system 
against DNS and concluded that it is suffient to predict 
both the linear instability and the early nonlinear evolu- 
tion of this instability. We showed that DNS agrees very 
well with the linear predictions in all cases. In many 
cases, particularly for small M and in the stable con- 
figurations, the four-wave truncation predicts quite well 
the early and intermediate nonlinear evolution phases. 
The nonlinear evolution for small M's is characterized 
by dominant wave dynamics, whereas for large M's the 
nonlinear evolution leads to rolling up of the carrier wave 
vorticity into Karman-like vortex streets. Such hydrody- 
namic vortices behave very differently from waves, and 
it it precisely at the moment of roll-up that the full sys- 
tem's evolution strongly diverges from the prediction of 
the four-wave truncation. After the roll-up the full sys- 
tem enters into a saturated quasi-stable state which per- 
sists for a relatively long time but eventually decays due 
to presence of hyper-viscosity. On other hand, the corre- 
sponding four-wave system keeps going though an infinite 
sequence of nonlinear oscillations. If M is small and the 
roU-ups do not occur (or are delayed) the full system may 
follow its four-wave counterpart for much longer: its ini- 
tial growth can reverse and may exhibit the nonlinear 
oscillations associated with the 4MT. 

Finally, we would like to emphasize two physical effects 
that can be important for both plasma and GFD systems. 
For M > 1 we observe the formation of stable, narrow 
zonal jets, in agreement with earlier theoretical predic- 
tions of As we mentioned, these jets are more stable 
than one would expect based on the Rayleigh-Kuo crite- 
rion alone because their 2D structure consists of stable 
vortex streets. Such narrowjets represent very effective 
transport barriers which may be responsible for the LH 
transitions in tokamaks. Indeed, narrowness of the jet 
causes characteristic pedestal-like radial profiles of the 
particles and energy typically seen in tokamaks during 
H-mode, i.e. the narrow jet provides a thin "insulation" 
layer (called internal transport barrier) which keeps the 
particle and energy densities significantly higher inside. 
The second physical effect we would like to mention oc- 
curs at low nonlinearities M. This is the fact that the sys- 
tem tends to select the states with somewhat off-zonal jet 



structures. This tendency to favour the off-zonal struc- 
tures is seen already on the level of the linear analy- 
sis, where as we showed the most unstable modulation 
changes from zonal to off-zonal when the nonlinearity 
is reduced. Possibly, this mechanism can explain recent 
ocean observation of off-zonal jet striations [7] (note that 
the nonlinearity of the ocean Rossby waves in these sit- 
uations is likely to be rather low). 

In future, it would be interesting to study in more de- 
tail how the transport properties are affected by both 
zonal and off-zonal jets that arise in the strongly and 
weakly nonlinear cases respectively. In particular, it 
would be interesting to see how different are the trans- 
port barriers provided by coherent vortex streets from 
the barriers provided by random jets at the later stages. 
It would also be interesting to study situations where a 
broad spectrum of modulations is present initially and, 
in particular, to verify that the system selects the most 
unstable one. If the carrier wave spectrum is not narrow 
it would be of interest to study when the modulational 
instability wins over the inverse cascade mechanism. 
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